Size-dependent spinodal and miscibility gaps for intercalation in nano-particles 
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Using a recently-proposed mathematical model for intercalation dynamics in phase-separating 
materials [Singh, Ceder, Bazant, Electrochimica Acta 53, 7599 (2008)], we show that the spinodal 
and miscibility gaps generally shrink as the host particle size decreases to the nano-scale. Our work 
is motivated by recent experiments on the high-rate Li-ion battery material LiFeP04; this serves 
as the basis for our examples, but our analysis and conclusions apply to any intercalation material. 
We describe two general mechanisms for the suppression of phase separation in nano-particles: (i) a 
classical bulk effect, predicted by the Cahn-Hilliard equation, in which the diffuse phase boundary 
becomes confined by the particle geometry; and (ii) a novel surface effect, predicted by chemical- 
potential-dependent reaction kinetics, in which insertion/extraction reactions stabilize composition 
gradients near surfaces in equilibrium with the local environment. Composition-dependent surface 
energy and (especially) elastic strain can contribute to these effects but are not required to predict 
decreased spinodal and miscibility gaps at the nano-scale. 



Introduction. Intercalation phenomena occur in 
many chemical and biological systems, such as graphite 



intercalation compounds [ij 



DNA molecules 



solid- 



oxide fuel cell electrolytes [3|, and Li- ion battery elec- 
trodes [4J. The intercalation of a chemical species in a 
host compound involves the nonlinear coupling of surface 
insertion/extraction reaction kinetics with bulk transport 
phenomena. It can therefore occur by fundamentally dif- 
ferent mechanisms in nano-particles and molecules com- 
pared to macroscopic materials due to the large surface- 
to-volume ratio. Intercalation dynamics can also be fur- 
ther complicated by phase separation kinetics within the 
host material. This poses a challenge for theorists, since 
phase transformation models have mainly been developed 
for periodic or infinite systems in isolation Q , rather than 
nano-particles driven out of equilibrium by surface reac- 
tions. 

In this paper, we ask the basic question, "Is nano differ- 
ent?" , for intercalation phenomena in phase-separating 
materials. Our analysis is based on a general mathemat- 
ical model for intercalation dynamics recently proposed 
by Singh, Ceder, and Bazant (SCB) @. The SCB model 
is based on the classical Cahn-Hilliard equation f?\ with a 
novel boundary condition for insertion/extraction kinet- 
ics based on local chemical potential differences, includ- 
ing concentration-gradient contributions. For strongly 
anisotropic nano-crystals, the SCB model predicts a new 
mode of intercalation dynamics via reaction-limited non- 
linear waves that propagate along the active surface, fill- 
ing the host crystal layer by layer. Here, we apply the 
model to the thermodynamics of nano-particle intercala- 
tion and analyze the size dependence of the miscibility 
gap (metastable uniform compositions) and the spinodal 
region (linearly unstable uniform compositions) of the 
phase diagram. 



Our work is motivated by Li-ion battery technology, 
which increasingly involves phase-separating nanoparti- 
cles in reversible electrodes. The best known example is 
LiFeP04, a promising high-rate cathode material @ that 
exhibits strong bulk phase separation . Experi- 

ments have shown that using very fine nano-particles (< 
100 nm) can improve power density and (with sur- 

face modifications) achieve "ultrafast" discharging of a 
significant portion of the theoretical capacity jl 3 | . Exper- 
iments also provide compelling evidence [id. Il4 [TBI , [lij 
for the layer-by-layer intercalation waves (or "domino 
cascade" [l3]) predicted by the SCB theory 0, [l3|, in 
contrast to traditional assumption of diffusion limitation 
in battery modeling [l^, [l^ . 

There is also experimental evidence that the equilib- 
rium thermodynamics of LiFeP04 is different in nano- 
particles. Recently, Meethong et al. have observed that, 
as the crystal size decreases, the miscibility gap be- 
tween the lithium-rich and lithium-poor phases in the 
material shrinks significantly [20] (i.e. the tendency for 
phase separation is reduced). A suggested explanation 
is that smaller particles experience relatively larger sur- 
face effects, which has been supported by calculations 
with an elaborate phase-field model [2l[, although with- 
out accounting for surface reaction kinetics. The size- 
dependency of the miscibility gap is fairly strong in 
LiFeP04, and the importance of the surface energy has 
been demonstrated. However, it has also been seen ex- 
perimentally that carbon coating can reduce the surface 
effects and prevent the surface-induced reduction of the 
miscibility gap [2^ . 

We will show that the SCB model suffices to predict 
that the spinodal and miscibility gaps both decrease as 
the particle size decreases. The analysis reveals two fun- 
damental mechanisms: (i) nano-confinement of the inter- 
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phase boundary, and (ii) stabilization of the concentra- 
tion gradients near the surface by insertion/extraction re- 
actions. These effects are independent of surface energy 
models, and indeed are valid for any phase-separating 
intercalation system. 

Model. We employ the general SCB model for inter- 
calation dynamics — based on the Cahn-Hilliard equation 
with chemical-potential-dependent surface reactions — 
without any simplifying assumptions In particular, 
we do not specialize to surface-reaction-limited or bulk- 
transport-limited regimes or perform any depth averag- 
ing for strongly anisotropic crystals 0, • Our field of 
interest is c(x, t) , the local concentration of the interca- 
lated diffusing species (e.g., Li in LiFeP04). Let p be the 
density of intercalation sites per unit volume in the sys- 
tem (e.g., occupied by Li ions or vacancies), assumed to 
be constant and independent of position and local con- 
centration. We take c to be normalized by p, so it is 
non-dimensional and only takes values between and 1 
(e.g., in the local compound LicFeP04). 

We assume that the free energy of mixing in our model 
system is well-approximated by the Cahn-Hilliard func- 
tional SSEl 



ghom(c) + -(Vc)^K(Vc) 



pdV . (1) 



The function ghom (c) is the free energy per molecule of a 
homogeneous system of uniform concentration c, which is 
non-convex in systems exhibiting phase separation. The 
gradient penalty tensor K is assumed to be a constant 
independent of x and c. Then the diffusional chemical po- 
tential (in energy per molecule) is the variational deriva- 
tive of Gniix, 



^(X,i) ^ ^ghom(c) _ ^ 

oc 



(2) 



The mass flux (in molecules per unit area per unit time) 
is given by the linear constitutive relation 24 1 



J(x,t) = -pcMVp , 



(3) 



where M is a mobility tensor (denoted B by SCB 
Finally, the dynamics are governed by the mass conser- 
vation equation 



djpc) 
dt 



V-J = 



(4) 



For illustration purposes, we employ the regular solu- 
tion model for the homogeneous free energy [251] : 

5hom(c) = ac(l - c) + ksT [clogc + (1 - c) log(l - c)] . 

The two terms give the enthalpy and entropy of mixing, 
respectively. When numerical values are needed, we will 
use a/fcsT = 5, which is in rough agreement at room 



temperature with measurements on LiFeP04 [2611. Of 
course, other models are possible, but for the intercala- 
tion of a single species in a crystal with bounded com- 
positions < c < 1 , the homogeneous chemical potential 
(c) = Shorn (c) must diverge in the limits c ^ 0+ 
and c — > 1~ due to entropic contributions from particles 
and vacancies. (This constraint is violated, for example, 
by the quartic ghom(c) from Landau's theory of phase 
transitions, suggested in a recent paper on LiFeP04 [l6l | 
following SCB.) 

Note that K/fc^T has units of length-squared. Since it 
is assumed that K is positive-definite, we may denote its 
eigenvalues by fcsTAf for real, positive lengths A^. In par- 
ticular, when K is diagonal, we define Xi = Ka/ksT . 
When the system is phase-separated into high-c and low- 
c regions, these \i are the length scales for the inter- 
phasial widths in the different eigendirections 0, [23|] . In 
LiFeP04, experimental evidence [lO|] suggests that one of 
these widths is about 4nm (though the in the other 
two directions might be large — comparable to the par- 
ticle size — as phase-separation in these directions is not 
believed to occur). These are therefore the natural length 
scales for measuring the size of phase-separating nano- 
crystals. 

Our system of equations is closed by the following 
boundary conditions on the surface of the nano-particle: 



n ■ (KVc) 
n J 







(6) 

-PsR (7) 

where fi is an outward unit normal vector. Equation [S] 
is the so-called variational boundary condition, which is 
natural for systems without surface energies or surface 
diffusion and follows from continuity of the chemical po- 
tential at the surface. Equation[7]is a general flux condi- 
tion enforcing mass conservation, where ps is the surface 
density of intercalation sites, and R is the net local rate of 
intercalant influx (insertion) across the boundary. In the 
classical Cahn-Hilliard (CH) model, no mass flux across 
the boundary is allowed, and thus i? = 0. For intercala- 
tion systems [l3| , we allow for a non-zero reaction rate 
R depending on the local values of c and p and refer to 
this general set of equations as the Cahn-Hilliard- iwit/i- 
reactions (CHR) system. 

For the current work (and indeed, for many of the con- 
clusions reached by SCB @), the particular form of R 
is unimportant. According to statistical transition-state 
theory in a concentrated solution, the net insertion rate is 
given by the difference of insertion and extraction rates, 
each having Arrhenius dependence on an (excess) chem- 
ical potential barrier. In order to satisfy de Bonder's 
equation 25[, it must have the general form 



R — Ru 



1 — exp 



P- Pe 



(8) 



where i?ins is the rate for the insertion reaction. For ther- 
modynamic consistency, this p must be the same as the 
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(a)Mixing energies of steady-state, phase-separated solutions. (b)Width of the miscibility gap as a function of crystal size. 



FIG. 1: The free energies in the first plot are given per intercalation site so that they are comparable across different crystal 
sizes. The dotted line indicates the free energy of mixing per site for a uniform system of the given concentration. 



diffusional chemical potential used in the bulk equations, 
and /ie is the external chemical potential of the inter- 
calants in a reservoir phase outside of the particle (e.g., 
Li+ in the electrolyte and e~ in the metallic current col- 
lector of a Li-ion battery electrode); note that we are 
again assuming that the particle surface is energetically 
identical to the bulk. If the reaction rates were controlled 
by electrostatic potential differences, for example, then 
i?ins could include transfer coefficients and the interfacial 
voltage drop, and the familiar Butler- Volmer model for 
charge-transfer reactions 27[ would be recovered in the 
limit of a dilute solution. Following SCB Q, we do not 
make any dilute solution approximation and keep the full 
CH expression for /i ^ — including the second derivative 
term — while assuming a uniform external environment 
at constant /ig. Although different models are pos sible 
for the chemical potential of the transition state [28[ , we 
make the simple approximation of a constant insertion 
rate i?ins, consistent with particles impinging on the sur- 
face at constant frequency from the external reservoir. In 
that case, the composition dependence of R enters only 
via the extraction rate. 

The CH Miscibility Gap. Outside of the spinodal 
range, systems with uniform concentration fields are lin- 
early stable. However, if there exists a phase-separated 
solution with the same overall amount of our material 
but with a lower free energy, then the uniform system will 
only be metastable. We will demonstrate that the misci- 
bility range — the set of overall concentrations for which 
phase separation is energetically favorable — shrinks as 
the particle size decreases. 

Unlike the spinodal, the miscibility gap cannot be stud- 
ied analytically. Instead, we must solve our original set of 
equations (QHll) numerically, looking for phase-separated 
systems with lower free energies than the uniform system 



with the same overall concentration. We focus only on 
1-dimensional systems, or equivalently 3-dimensional sys- 
tems whose phase boundary is perpendicular to one of the 
eigendirections and whose concentration field is uniform 
in the other two directions. Note that there is experi- 
mental [ifll and theoretical 15 1 evidence that this is an 
accurate picture for the concentration field in LiFeP04. 
We will henceforth drop the subscripts on A, and call L 
the length of the system. 

We begin by fixing a single crystal size. For each value 
of the average concentration, we choose a corresponding 
initial condition, and we solve the Cahn-Hilliard equa- 
tion (using a semi-discrete finite volume method and the 
no- flux boundary condition). The system is stepped for- 
ward in time until the free energy reaches a minimum. 
The resulting free energies of mixing for three different 
crystal sizes and a range of average concentrations are 



plotted in Fig. 1(a) Note that the curves do not extend 



across the entire x-axis. This is because, for sufficiently 
extreme average concentrations, no initial conditions can 
be found which lead to a phase-separated steady state. 
This suggests that such states do not exist, or that if 
they do exist they are not local minimizers of Gmix- The 
phase-separated energy curves do extend slightly past the 
uniform curve, allowing us to estimate the endpoints of 
the miscibility gap. The results suggest that the misci- 
bility gap shrinks as the crystal size decreases. 

In order to validate this hypothesis, we performed a 
more exhaustive search for phase-separated, steady-state 
solutions near the apparent miscibility endpoints. This 
was done using the shooting method for boundary value 
problems to compute concentration fields satisfying ^ 
with iL = constant. The resulting field that extremized 
the average concentration while still having a smaller free 
energy of mixing than the corresponding constant field 
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was considered to be the boundary of the miscibihty re- 
gion. The calculated miscibility gap widths over a range 



of crystal sizes are plotted in Fig. 1(b) they clearly sup- 
port a shrinking miscibility gap. 

There is a simple physical explanation for this behav- 
ior. As discussed above, the interphase region will nor- 
mally have a width on the order of A. The average con- 
centration can only be close to or 1 if this interphase 
region is close to a system boundary. At this point, the 
average concentration can only become more extreme if 
the interphase region is compressed or disappears. If it 
disappears, then we are left with a uniform system, and 
the average concentration is outside of the miscibility 
gap. The other alternative, though, is expensive ener- 
getically due to the gradient penalty term in ([T]). Thus 
low-energy, phase-separated systems are limited geomet- 
rically to those concentrations in which the interphase 
region is (relatively) uncompressed between the crystal 
boundaries. As the crystal size decreases, the limits im- 
posed on the average concentration by the incompress- 
ibility of the interphase region becomes more and more 
severe, and thus the miscibility gap must shrink. 

The CHR Spinodal Gap. The spinodal gap is 
the set of concentrations for which an initially-uniform 
system will spontaneously decompose through the expo- 
nential growth of infinitesimal fluctuations. Thus, per- 
turbation theory is the relevant mathematical tool, and 
we look for solutions to the CHR system of the form 

c(x,t) = Co + eci(x,t) 

where cq is a constant and e is a small parameter. If cq 
is truly a static solution to the CHR equations, then by 
([5]), fj,e must equal /io = 95hom(co)/9c at all points on 
the boundary of V. The first-order system derived by 
linearizing about cq is then 



Mi(x, t) 
Ji(x, t) 

d{pci) 
dt 



= -V- Ji 



ci - V • (KVci) 



with the boundary conditions 
n- (KVci) = 

PsRu 



n • Ji = 



(9a) 
(9b) 

(9c) 

(9d) 
(9e) 



This is a fourth-order, linear system with constant coef- 
ficients. Note that the exact same set of equations would 
result even had we taken K and M to be functions of c; 
the tensors above would only need to be replaced by the 
(still constant) values K(co) and M(co). 

If we have an infinite system with no boundaries, then 
the Fourier ansatz e^^'^e^*^ solves the above system if and 




FIG. 2: Width of the spinodal gap as a function of crystal 
size. The dotted line indicates the width of the spinodal re- 
gion for an infinite system. The other three curves are given 
for different values of the non-dimensionalized reaction rate 
constant 7^ = (ps/ p\){R\ns/ {D /\^)), where D = MksT is 
the difi'usion constant in the bulk. 



only if it satisfies the dispersion relation 

5 = -co(k^Mk)(^^^^|^-Hk^Kk) . (10) 

Since M and K must be positive-semidefinite, s will be 
non-positive whenever 9^(7hom(co)/9c^ > 0. However, if 
ghom{co) / dc^ < 0, then the cq will be unstable to long- 
wavelength perturbations. In particular, for the regular 
solution model ([5]), the criterion for linear stability be- 
comes 



+ 



1 



ksT co(l-co) 



> 



Thus a high enthalpy of mixing will promote instability 
of uniform systems with moderate concentrations. 

If instead the system geometry is finite, then the 
boundary conditions will constrain the set of allowable 
wave vectors k. We again focus on one-dimensional sys- 
tems for simplicity. Then if the system occupies the line 
segment from to L, the general solution of the per- 
turbed equations for the CH system (i?ins = 0) is a sum 
of terms of the form 



/ \ A Z'^"' 

ci(x,t) — Acqs y—- 



for any integer n. The dispersion relation pH]) still holds, 
but the wave number must equal m:/L for integer val- 
ues of n in order to satisfy the boundary conditions. In 
other words, we can no longer perturb the system with 
arbitrarily-long wavelengths. The stability criterion is 

For the regular solution model ([S]), the criterion for 
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linear stability becomes 



ksT co(l-co) 



The spinodal region is defined as the range (a, 1 — a) 
of unstable cq values. It is easily verified that a is a 
decreasing function of L, i.e. that the spinodal range is 
more narrow for smaller crystals. Moreover, for suffi- 
ciently small values of X/L, the above inequality is sat- 
isfied for all values of cq, in which case there is no spin- 
odal region at all. These facts are demonstrated in [21 
These results date back to Cahn's 1961 paper [3] and are 
known in the phase-field community. However, it seems 
that their relevance for nano-particle composites — as in 
Li-ion batteries — has not yet been appreciated. 

Moving beyond classical bulk models, we will now 
show that non-zero boundary reactions can further re- 
duce the spinodal gap width. Even the linear perturbed 
system of equations is no longer analytically tractable 
when i?ins 7^ 0, and in particular, the wave numbers are 
no longer simply mr/L. According to the dispersion rela- 
tion ()10p . every s is associated with four wave numbers, 
and in general it takes a linear combination of all four 
such functions to satisfy the boundary conditions. For 
any given L, cq, and s, we may compute the four cor- 
responding wave numbers kj, and then look for a set of 
coefficients Aj such that ^^e^'^^^e^* solves the per- 

turbed PDE and boundary conditions Because the 
system is linear and homogeneous, this can be reduced 
to finding a solution to some matrix equation BA = 0, 
which has solutions if and only if the determinant of the 
matrix B is 0. 

Therefore, for any given system size L and reaction rate 
constant i?ins, we must numerically solve for the range of 
concentrations cq that admit solutions to the perturbed 
equations for at least one positive value of s. Results of 
such computations are shown in [2l Notice that increas- 
ing the reaction rate constant reduces the spinodal gap. 
Moreover, it was found numerically that increasing i?ins 
tends to reduce the growth rate constant s. 

These effects cannot be explained solely in terms of 
chemical potential perturbations near the boundary. In- 
stead, we must examine the nature of the allowable per- 
turbations for different reaction rates. For large values of 
i?ins, any non-zero ni at the boundaries causes large per- 
turbations in the reaction fluxes by Eq. [9el In order for 
Ji to be diflerentiable near the boundaries, we must also 
have large bulk fluxes nearby. In general, this would re- 
quire large concentration gradients, or equivalently short- 
wavelength perturbations. But, as is clear from the dis- 
persion relation , it is precisely the s/iori- wavelength 
perturbations which are rendered stable by the gradient 
penalty term (see Ref. for an interesting discussion 
of this point). 

More mathematically, suppose fii is non-zero at a 
boundary. Then by ([5e)) . Ji must be non-zero there, 



which by (j9bp implies that /xi must have a non-zero gra- 
dient. Combining these two terms with our ansatz for ci 
yields the requirement 



Co 



{ksTM) k 



We therefore see that when fii is non-zero at a boundary, 
the wave number scales linearly with the reaction rate 
constant. Again, large i?ins would require large k, which 
are increasingly stable. 

As the reaction rate increases, then, unstable perturba- 
tions satisfying ([9]) must have fii and V/ii close to near 
the boundaries. However, this requires long-wavelength 
perturbations, and we have already shown that these will 
become increasingly stable as the crystal size shrinks. 
Thus fast reaction rates will tend to stabilize small nano- 
particles. 

Note, however, that for systems larger than about 2.5A, 
the spinodal gap does not disappear even for infinitely 
fast reactions. This implies that there must exist in- 
finitesimal perturbations to a uniform system which lead 
to phase separation without ever changing the difFusional 
chemical potential at the boundaries of the system. This 
has been verified numerically by solving the full CHR sys- 
tem ([SHE]) in the limit _Rins oo. However, by limiting 
the spinodal decomposition to only occur via this small 
class of perturbations, higher reaction rates reduce the 
decomposition growth rate and the spinodal gap width. 

Though we have used a specific mathematical model to 
derive these results, the conclusions are generally valid. 
Regardless of the bulk model, a bounded system will only 
allow a discrete spectrum for its first-order perturbations. 
The smallest admissible wave numbers will scale like 1/ L, 
and the system will suffer linear instability for more nar- 
row ranges of concentrations as the system size shrinks. 
Moreover, fast reaction rates at the boundaries require 
short wavelength perturbations, and such perturbations 
are energetically unfavorable when there is a diffuse in- 
terface between phases. 

Other Effects. There are at least two important 
effects which we have excluded from our analysis. First, 
we have intentionally ignored surface energies in order to 
demonstrate that purely bulk effects and reaction rates 
can cause shrinking spinodal and miscibility gaps. How- 
ever, surface energies could easily be accommodated. For 
example, if the free energy of the system were to include 
a concentration-dependent surface tension between the 
particle and its environment, 

Gmix = G,„ix,bulk + / 7(c) dA , 
JdV 

then the variational boundary condition ^ would need 
to be replaced by 



n • (KVc) 



1 dj 
p dc 



. 
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This would change the analysis, but would not signifi- 
cantly affect the conclusions. 

Perhaps a more serious omission for LiFcP04 in par- 
ticular is elastic stress in the crystal due to lattice mis- 
matches. However, it has been demonstrated [s^l that 
these effects can frequently be accommodated by simply 
decreasing the enthalpy-of-mixing parameter a. Given 
our results above, elastic stress would therefore enhance 
the shrinking spinodal and miscibility gaps. 

Conclusion. We have shown that intercalation phe- 
nomena in phase-separating materials can be strongly 
dependent on nano-particle size, even in the absence of 
contributions from surface energies and elastic strain. In 
particular, the miscibility gap and spinodal gap both de- 
crease (and eventually disappear) as the particle size is 
decreased to the scale of the diffuse interphase thickness. 
Geometrical confinement enhances the relative cost of 
bulk composition gradients, and insertion/extraction re- 
actions tend to stabilize the concentration gradients near 
the surface. These conclusions have relevance for high- 
rate Li-ion battery materials such as LiFeP04, but are in 
no way restricted to this class of materials. 
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